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ABSTRACT 

We perform a linear stability analysis for corrugations of a Newtonian shock, 
with particle pressure included, for an arbitrary diffusion coefficient. We study 
first the dispersion relation for homogeneous media, showing that, besides the 
conventional pressure waves and entropy /vorticity disturbances, two new pertur- 
bation modes exist, dominated by the particles' pressure and damped by diffu- 
sion. We show that, due to particle diffusion into the upstream region, the fluid 
will be perturbed also upstream: we treat these perturbation in the short wave- 
length (WKBJ) regime. We then show how to construct a corrugational mode 
for the shock itself, one, that is, where the shock executes free oscillations (pos- 
sibly damped or growing) and sheds perturbations away from itself: this global 
mode requires the new modes. Then, using the perturbed Rankine-Hugoniot 
conditions, we show that this leads to the determination of the corrugational 
eigenfrequency. We solve numeri cally the equation s for the eigenfrequency in the 



WKBJ regime for the models of lAmato and Blasil . showing that they are stable. 



We then discuss the differences between our treatment and previous work. 
Subject headings: shock waves - cosmic rays 



1. Introduction 

One of the major uncertainties still surrounding particle acceleration around shocks is 
the level at which this saturates: i.e., the total energy content in non-thermal particles, 
per unit volume, TZ, in units of the incoming fluid kinetic energy density. This saturation 
level plays a very important role, for instance, in discussions of the origin of cosmic rays as 
observed at Earth: it has to be rather large (> 0.1) to allow SNe to provide the observed 
flux. Also, discussions on the origin of UHECRs are influenced by similar considerations: 
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the idea that UHECRs are accelerated by GRBs (jVietril Il995l ; IWaxmanl Il995l ) draws some 
measure of support by the coincidence between the energy rel ease rate of GRBs in 7-ray 



photo ns and that required to account for observations at Earth (I Waxmanl |2002| ; IVietri et al 
2003h . 



It is certainly possible that this level of saturation is somehow connected with the very 
uncertain particle injection mechanism, and that saturation occurs, in many astrophysically 
important situations, at very low levels. Yet, given the arguments concerning SNe and GRBs 
mentioned above, it appears that at least occasionally saturation must occur at rather large 
levels. This is especially so in the case of Galactic cosmic rays: if in fact we were to discard 
SNe as sources of cosmic rays, all other possible classes of sources, being both less numerous 
and less energetic, would force us to require saturation at even larger levels^. 

A distinct possibility is that the saturation level is determined not by the injection mech- 
anisms, but by modifications which the particles' pressure induces in the shock properties. 
After all, in the test particle limit the particles' spectrum is ultraviolet divergent (though 
marginally so, of course), and one may hope that, since the particles' back-reaction on the 
fluid makes acceleration less efficient (by reducing the velocity jump around the gas sub- 
shock V_a_convergent spectrum will be obtained, with a definite value for the parameter 1Z. 
After iMalkovl 's seminal papers jMalkov! Il997l : iMalkov et all I2OO0I : iMalkov and O'C. Drurv 



200ll ). ful ly self-consistent solutio ns with particles' pressure properly included have been ob- 
tained by lAmato and Blasil (120051 ). In these models, it appears that particles can carry away 
an arbitrarily large fraction of the incoming energy flux, for sufficiently large Mach number 
at upstream infinity. Even more intriguing is the fact that these solutions have particles' 
spectra which are more, not less, ultraviolet divergent than those in the test particle limit. 
In fact, this divergence is arrested only by arb itrarily limiting the largest individual energy 
to some fiducial value (lAmato and Blasil 120051 ). 



These solutions then seem to beckon for a stability analysis, in the hope that they are 
found to be unstable once 1Z exceeds some critical value. The instability we envision is of 
course the corrugational instability for shocks, whereby ripples on the shock surface become 
of larger and larger amplitud^E If this instability were to exist, we might easily imagine that 
the shock is substituted by a region of strong fluid turbulence, where particles moving a few 
Larmor radii perceive only a small velocity difference at the two ends of their free wandering, 



2 The situation is much less clear for UHECRs, because many of the proposed classes of sources are purely 
hypothetical, so that there are no constraints on their properties. 

3 In all of this paper, the term corrugation instability is taken to include also the spontaneous emission of 
sound waves by the shock. 



- 3- 



and thus acceleration to high energies is made somewhat less efficient. 



However, it is well-known (ILandau and Lifshita 119871 ) that shocks in polytropic fluids 



are stable against this kind of perturbations, and against the spontaneous emission of sound 
waves as well. The hope for the existence of an instability is based on a rather more subtle 
argument. When the shock surface flaps, it sheds in the downstream region pressure waves 
as well as entropy and vorticity perturbations. We will show later that the last two do not 
couple to particle perturbations, but pressure waves do, thus generating small perturbations 
in the particles' distribution function, Sf. These particles will however return upstream by 
means of diffusion, and here generate, by their sheer perturbed pressure, some more fluid 
perturbations. Thus, not all energy shed by the shock is lost forever toward upstream infinity, 
as is the case in the purely hydrodynamical case, but some fraction of it returns to the shock 
to generate more flapping. In fact, since pressure waves are strongly damped by diffusion, 
and nearly all particles return to the upstream section since the shock is Newtonian, we may 
conjecture that most energy shed by the shock flapping makes it back to feed more flapping, 
even though account must be taken of diffusion and phase mixing. Is it possible that this 
sets up some strong reinforcement, making the whole system unstable? This is the question 
we address in this paper. 

This question has of course been studied before, under somewhat different conditions. 
The diffusion coefficient of non-thermal particles by the fluid (= D) has been taken (by 
other authors) as independent of the particle impulse: this is conventionally referred to as 
the two-fluid approximation, because the Boltzmann equation for the non-thermal parti- 
cles can then easily be recast into an equation for their pressure, thusly erasing all refer- 
ences to the underlying distri bution function. Besides making the two-fluid approximation, 
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(19 98) a l so neg lected diffusion altogether. A more complex analysis 



has been presented by lToptyginl (119991 ). who included energy transport and particle injection 
at the shock, but still in the two-fluid approximation. In the following, we shall make no such 
approximations: we will consider a finite diffusion coefficient D, and it will be allowed to 
be an arbitrary function of both p, the particle momentum, and of p, the local fluid density, 
so as to at least mimic the increase in magnetic field strength due to flux freezing. The 
arbitrary nature of the dependence of D on p automatically prevents the use of the two-fluid 
approximation, and forces us to use the full Boltzmann- Skilling equation. 

A word of caution is in order about some assumptions. We will neglect all energy in the 
form of magnetic field and Alfven waves; of course this is necessary because the zero-th order 
solutions of Blasi and collaborators do not include these effects, but in our case this neglect 
requires one extra assumption, i. e. that the time scales for energy to accumulate into any of 
these energy sinks, T in , and to flow out of each of them, T out , be ordered like this: T out < T in . 
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If this inequality were severely violated, the magnetic field or Alfven waves might acquire a 
significant fraction of the total energy, despite their negligibility in the zero-th order solution, 
and make our treatment completely irrelevant. 

In the foll owing, we shall follow closel y the treatment of the shock corrugational insta- 
bility given by lLandau and Lifshita (119871 ). In particular, we shall consider a steady-state 
shock in its own frame, located at z = 0, with fluid coming from the left , and exiting from 
the right, so that all speeds are > 0. Exactly like lLandau and Lifshitzl . we shall consider 
perturbations generated by the shock flapping, so that there can be no incoming waves, from 
either upstream or downstream infinity. In Section [2J we shall consider perturbations in a 
homogeneous medium; this is perhaps a tad boring, but it contains a number of new results 
which are of the utmost importance later on. In Section [3] we briefly discuss perturbations 
upstream, i.e. where the flow is inhomogeneous; we show here that we can easily obtain 
the perturbations in the WKBJ limit k y — > oo, which restricts our analysis to the regime 
A < L, where A is the perturbation wavelength perpendicular to the shock, and L is the 
typical size of the region of inhomogeneity upstream. In Section H] we discuss the boundary 
conditions on the perturbed particle distribution function, and derive how it relates to the 
amplitude of the modes to which particle perturbations couple. We present in Section \5\ the 
perturbed Rankine-Hugoniot conditions and in Section [6] what fixes the global corrugation 
mode eigenfrequency. In Section [7] we present our nu merical computations for the stability 
of the exact solutions of the zero-th order problem by lAmato and Blasil (120051 ) . In Section 
[HI we will compare our results with other works in the literature, and briefly summarize our 
work. 



2. Perturbations in a homogeneous medium 

We give here, for future reference, our basic equations. They are the conventional 
hydrodynamic equations: 

^ + V " (pv) = (1) 

J^?+tf. V ir=_I v (p + p c ) ( 2 ) 

3s 

— + v. V s = (3) 

which contain a term for the momentum exchange between the fluid and the non-thermal 
particles represented by the gradient of the particle pressure P c , plus the usual Boltzmann 



- 5 - 



equation in ISkillingl (119751 ): 

^ = V • (D V /) - v ■ V f + \ V -vp^ ■ (4) 

We assume D = D(p, p) to be a given function of p and p. 

We consider small-amplitude perturbations around a homogeneous solution where the 
particles are supposed to exert a non-negligible pressure. First, we consider entropy pertur- 
bations. Perturbations can be taken in the form 



5s oc exp ynut — %k ■ f) (5) 

so to obtain, from the equation of entropy conservation, eq. [3J 

{uj — uk x )5s = — > uj — uk x = , (6) 

where u is the unperturbed fluid velocity in the x-direction. Perturbation of the mass 
conservation equation yields 

{uj — uk x )Sp — pk-Su = 0^k-Su = 0. (7) 

Perturbation of the momentum equation yields 

(uk x - uj)Su + -(SP + 5P C ) = -> SP = -6P C . (8) 
P 

Now, the perturbation of the Boltzmann equation yields 

iu5f = -k 2 D5f + mkjf - -k ■ Sup^- (9) 

— * 

which implies 5f = because, for entropy perturbations, u = uk x and k ■ Su = 0, as deduced 
above. Furthermore, since 5f = 0, necessarily 5P C = 0, and thus SP = 0. It follows that 
entropy perturbations and particle perturbations are completely decoupled: in fact, entropy 
perturbations are just advected by the zero-th order flow. 

In summary, entropy perturbations have the following characteristics, where we define, 
for future convenience, V = 1/p: 

uj — uk x = 

SP = SP C = 8f = SI = 

SV 7 — I m6s , . 

(10) 



V 7 k 



B 
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where the last equation applies to ideal fluids: m is the average particle mass and ks Boltz- 
mann's constant. 

We turn now to isentropic perturbations, which automatically satisfy eq. [3J Mass 
conservation implies 

{ijj — uk x )5p = pk ■ Su (11) 
while momentum conservation implies 

(u - uk x )5u = -(SP + 5P C ) . (12) 
P 

Let us begin by assuming that the left- hand- side of eqJTTl vanishes; the same is then true for 
\/ -Su. Multiplying eq. IT2"lby kA, we see that \/ ASu = 0, unless u = uk x . If the curl vanishes, 
then so does Su, because any vector with vanishing divergence and curl is a constant, which 
can always be set to zero by a suitable choice of reference system. So we must have uj = uk x 
for perturbations with non-zero vorticity; we thus find that SP + SP C = 0, and then, again 
from eq. [H that Sf = 0, implying SP = SP C = 0: vorticity perturbations do not couple to 
particles either. We have for them: 

S s = Sp = SP = Sf = SP C = 
uj — uk x = 
k-Sv = 0. (13) 

We consider now those perturbations where Sf does not vanish. We can show that here 
too there are two distinct classes of modes: in the first one, Sf is not coupled to the fluid 
quantities, while in the second one it is through the term k ■ Su ^ 0. The first class of modes, 
which we call d-mode, cannot be obtained directly from eq. for reasons to be made clear 
shortly. We use instead the correct form 

f = v.(^/)-«f , (u) 

where we dropped the term oc k ■ Su in keeping with our desire to find a solution totally 
uncoupled from the fluid. The solution of this equation is known from elementary courses: 
if (j) (x,y,p) is the initial distribution function at time t = (possibly dependent upon p), 
the solution at later times for u = is 

/f 1 (x~x ) 2 (y-yp) 2 

dxo / 4/o M x o,Vo,p)^-^e *m e *m , (15) 

and the solution for u ^ is just Sf(x — ut,y,p). At the same time, we must make sure 
that this solution does not ruffle the fluid: this obviously requires SP C = at all times. Now 
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integrating eq. [IHover Anp 3 vdp/3 we find 

85P C 4tt d 2 f 3 0dP c 

-W = Yd^J DWfpvt + u-^-, (16) 

which clearly shows that, in order to have 5P C = everywhere at all times we need 8P C = 
J D5fp 3 v dp = everywhere at the initial time. If we now Fourier-expand the initial condition 
(f) Q with respect to x, y, the above conditions become 

<f) (x,y,p) = J dk x J dk y g(p, k)e~ tkxX e~ tkyV , J g(p,k)p 3 v dp = J g(p,k)D(p)p 3 v dp = . 

(17) 

This completes the derivation of this purely damped d-mode, which will not perturb the 
fluid. Though it may look at this stage like a mathematical oddity, this mode plays a key 
role in the matching of boundary conditions between the upstream and downstream regions. 
It is also worth remarking why it cannot be derived from its Fourier- analyzed counterpart: 
the solution in question contains a term oc e~ ly/ *, which does not have a Fourier transform 
with respect to t. 

When 5P C ^ 0, velocity perturbations Sv ^ are induced in the fluid by the non- 
vanishing particles' pressure gradient, and eq. [9] can be solved for Sf as 

Sppdf i(uj-uk x ) 
1 P 3dpi{u- uk x ) + k 2 D ' 1 } 

If we integrate this over 47rvp 3 dp/3, we find 

5p f PM ^ 4 df -i{u-uk x ) 

oPc = — -rrdpvp - — , 19 

P J Pm 9 dp i(u - uk x ) + k 2 D 

where we assumed the non-thermal particles to have a minimum (p m ) and a maximum (pm) 
momentum. It is convenient to introduce a weighted diffusion coefficient D defined as follows: 

l r^vir^ 



l-zD(z) JPMfcfrvpiZL 



which shows D to be a function of z only. 



The above can be simplified a bit by integration by parts. If the integral were to extend 
from to +oo, we would have 



'o 
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where P c is the particle pressure in the unperturbed flow, and 7 C , the effective particle 
polytropic index, satisfies 4/3 < 7 C < 5/3. Since however the integral extends only from p m 
to pm, we define 

i-pm fit 

w|^:p.. (22) 

' J pm 1 



Thus eq. [19] can now be rewritten as: 



sp c = l' c P c --. Tk^T^ • (23) 



We remark that both 7^ and .D are obtained by weighing the zero-th order solution, so 
that they can be immediately computed as soon as this solution is available. 

We can now eliminate k ■ Su between eqs. [11] and [IS, and then use SP = cj,5p (with c s the 
sound speed because we are considering isentropic perturbations) and the equation above to 
obtain 

^ = l + ~ ^-—^ (24) 



-yP j lk ' 2 D( lk2 

uj—ukx ^ u—uk x 

where we have called 



n = u ~ uK (25) 

the comoving eigenfrequency, in suitably scaled units. 

Eq. [2H together with the definition of D (eq. I20l) . is the sought-after dispersion relation 
for the coupled small perturbations in a homogeneous medium. 



2.1. Properties of coupled modes 



The case where D is independent of p has been derived before (jPtuskinl Il98ll ) and 
coi ncides with the above equation. Surprisingly, the existence of this mode was not noticed 
by iToptyginl (119991 ). even though a careful treatment of his equations yields precisely the 
same dispersion relation as above; this neglect of this mode has important consequences to 
be discussed later on. 



It is best to begin our discussion with the case when D, and thus D, is a constant, 
independent of p. The eq. [2H then reduces to 



ikD 



){Q 2 -1) = Q 



7P 



(26) 
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which is a simple polynomial equation of the third order. In this case two modes reduce to 
pressure waves, as is most easily seen in the test-particle regime P c = 0. There is however 
a third solution which is only slightly more mysterious: in the test particle regime these 
modes represent a local over-density of particles dissipated by diffusion. When the test 
particle regime does not apply, a particle contribution to the sound speed is introduced by 
the term oc P c . This new mode (which we call the third mode) is the equivalent of the 
d-mode when however the conditions [17] are not satisfied: the basic idea is still the same, the 
particle overdensity is dissipated, but since the particle pressure does not vanish, the fluid 
is consequently ruffled. Notice also that there are two third modes, traveling in opposite 
directions. 

The situation is slightly more complex when D = D(p), because one must solve simulta- 
neously eqs. [23] and [201 We begin by considering the limit k — ■> 0. In this case, and assuming 
Q — > a constant, we easily find, to leading order in k, Q 2 = 1 + ■y' c P c /('-yP): as it must, the 
dispersion relation allows pure pressure waves, with the particle pressure providing a correc- 
tion to the (pure gas) sound speed, because for large perturbation wavelengths particles are 
entrained by the perturbation. We remark that, in this limit, pressure waves are supersonic, 
in the sense that they are faster than pressure wav es propaga t ing in a pure gas of the same 



thermodynamical state, a result already noticed by iToptygin] ( 119991 ). 



Assuming Q — > a constant, we lost a solution, so we search for the third mode solution 
in the form Q = ak+ higher order terms in k. We find: 

_ %kb(%k 2 i{yj-uk x )) i 

Now we use the definition z = ik 2 /{u — uk x ) = tk/(Qc s ) to simplify the above to 

1 - zD(z) = (28) 
which can be now used with eqs. [20] and [22] to obtain 



Pn 



where of course df /dp < 0. As a function of real z, the right-hand side above (where it exists) 
is easily seen to be a monotonically decreasing function of z, vanishing for z — > ±oo. In any 
realistic physical problem the integral must extend from a minimum (p m ) to a maximum 
momentum pm] since D(p) is expected to be a monotonically increasing function of p, we see 
that the integral above always exists for z < l/D(p M ) = 1/D M , and z > 1/D(p m ) = 1/D m , 
and it diverges exactly at z — and z = 1/D m . Thus the integral on the right-hand 
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side of the equation above spans the whole range from to — oo as z varies between — oo 
and +1/D(p m ), and the range +00 to as z varies between l/D m and +00. Somewhere in 
the range l/D m < z < +00 there is the one and only solution of the above equation. An 
illustration of the integral on the right-ha nd side of the previous equation is shown in Fig. 



[TJ for a specific distribution function from lAmato and Blasil (120051 ): the qualitative features 



of this plot are generic to all distribution functions we have tried. 

Furthermore, since z = ik 2 /(uj — uk x ), and the small perturbations were assumed to 
vary as e t( , u) - uk *) t ( the result that z > implies that all modes are damped by diffusion, as 
physical intuition obviously suggests. 

The discussion in the opposite limit, k — ► +00, is similar. If we assume Q — > a constant, 
we find the solution 

n 2 = 1 , (30) 

without the correction to the sound speed due to the presence of the particles' pressure: in 
the limit k — > +00 particles escape by free streaming, and do not contribute to the sound 
speed. 

Again, we lost a solution, so we now search for the third mode as Q = ak+ lower order 
terms in k. We obtain now: 

ikD(ik 2 1 (u — uk x )) 



n = — i — - — (31) 



which can be rewritten as 



zD(z) = 1 . (32) 

Comparing this with eq. [20], we see that the value of z we are searching for is the one that 
makes the integral on the numerator of the right-hand side of eq. [20] diverge. Following the 
previous discussion, we see that this is 

' (33) 



D(p m ) ' 

This is always positive, so that the solution is always damped by diffusion. This result has 
a simple physical explanation: when a small overdensity of particles is generated locally, the 
time-scale for damping of this overdensity is dictated by diffusion of the slowest particles. 

Again for illustrative purposes, the real and imaginary parts of sonic and third mode s 



are displayed in Fig. [2] for a specific distribution function from lAmato and Blasil (120051 ) . 



Again, the qualitative features are generic to all models we investigated. 

In summary, we have seen that the introduction of particles modifies the modes of a 
homogeneous medium by adding two new modes, one coupled and one uncoupled to the 
fluid, both strongly damped by diffusion. 
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Fig. 1. — The right-hand-side o f eq. I29l as a function of z on the real axis, for one numerical 
solution from Amato and Blasi ( 20051 ); here, p m = 0.003mc, pm = 10 5 mc, D{p) oc vp. 
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Fig. 2. — Real and imaginary parts for the three solutions of eq. EH as a function of 
perturbation wave-number k, for the same distribution function as in Fig.l and 7e-P c / (jP) ~ 
158.629 (upstream pressure at 0~). The pressure waves are seen to be damped with the 
same rate (and opposite oscillation frequencies), while the third mode is pure imaginary, 
corresponding to a pure damping. 
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For future reference, we give the expressions for all small quantities, in units of 8P C : the 
link between 5P C and 5p is obtained from Eqs. [20] and [22] an d the others follow easily. For 
later convenience, we use again V = 1/p rather than p: 

fop — ftp e (tuJt-ik y y-ik x x) 

SV 1 8P C ( ikD\ 5P C _ 5P C 

„ kc a SV ( j c P c V 1 \ k sv _sv 

\ s c s Q / 

p = ~ly , (34) 

where the definitions of f and q will come handy later on. In these equations, k x must be 
regarded as a known function of u and k y , specified by eq. 



2.2. Initial value problem 

As a preparation for later work, we discuss the initial value problem. This has some 
interest because perturbations in 5f belonging to the various modes are not mutually orthog- 
onal, and thus it may appear that initial conditions, especially when given only in terms of 
5f, cannot be decomposed into mutually independent modes. To fix ideas, let us consider a 
homogeneous zero-th order solution where all fluid quantities are unperturbed, but a small 
perturbation Sfo(x,y,p) at time t = is given: what are the amplitudes of the four modes 
that will be excited (two pressure waves, the third mode, and the d-mode)? First of course 
we Fourier- analyze 5fo, calling a(p) its amplitude. We must then have 

a(p)=A i 8f i + g(p), (35) 

where we have introduced some notation that will be useful in the following: a summation 
convention over i is understood, the A^s are the amplitudes of the pressure waves and third 
modes for i — 1,2, 4, respectively; g(p) is the amplitude of the d-mode, as in eq. [TTJ, and the 

Sfi's are 

1 df %{uJi — uk xi ) 

3 dp ifa - uk xi ) + kfD{p) ' 
Comparing this with eq. HH it becomes clear that the mode amplitudes Aj's are simply 
Spi/p. Here the quantities u>i, k X i, kf are supposed to be linked by the appropriate branch of 
the dispersion relation. For g(p) to represent a proper d-mode, we know that it must satisfy 
the two constraints in eq. [T71 thus we derive two conditions on the mode amplitudes: 

J p 3 va(p) dp = Ai J p 3 v5fi dp (37) 
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j D(p)p 3 va(p) dp = Ai J D(p)p i v5f i dp . (38) 

The last condition can be obtained by satisfying the requirement that the perturbations to 
both density and velocity vanish at the initial time. With reference to eq. [MJ we see that 
the two conditions, the vanishing of the density and of the velocity, imply 

J2A = 0, JV4< = 0, (39) 

i 

which are two more linear relations which, together with eqs. [37] and [381 determine the A^s. 
This simple example illustrates the importance of the d-mode, a kind of elephants' graveyard 
because, once particles join it, the perturbations they generate can only be dissipated away 
without ruffling the fluid ever again. 



3. Perturbations upstream 

It is useful to remark that, upstream, entropy perturbations are not allowed for arbi- 
trarily dishomogeneous flows, while the same is not true for vorticity perturbations. 

In fact, the perturbation of the entropy equation is: 

d5s .d5s n , An . 
_ + „(*)_ = (, (40) 

where we could neglect the term Sv x ds/dx because the fluid is isentropic in the unperturbed 
state. From the above, we see that entropy perturbations are advected by the background 
flow all the way from upstream infinity to the shock; we cannot however accept this, since in 
our problem all perturbations must have as a source the shock flapping. Perturbations riding 
all the way from upstream infinity belong to perturbations in the boundary conditions, and 
are thus irrelevant. Thus the perturbed fluid will be assumed adiabatic, from now on. 

The same argument does not apply to vorticity perturbations when the fluid is stratified, 
because they do couple to particle perturbations: in fact we easily obtain from eq. [2] that 
the equation for the vorticity, if = V A v, is: 

(l + HHH i7 v VpAV(p+iy - (41) 

In the absence of particles, this equation tells us that vorticity is exactly (i.e., not just to 
zero or first order) Lie-advected by the flow because, for adiabatic flows, Vp A VP = 0. But 
in the presence of particles and of spatial gradients it is easily seen that the particle pressure 
P c acts as source of vorticity perturbations. 
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In order to make progress, we compute the spatial dependence of the various modes 
upstream in the WKBJ approximation: i.e., in the limit k y — > oo, and in particular k^ 1 <C L, 
the size of the upstream precursor. In other words, we take all perturbed physical quantities 
to be of the form 

5X « Q x {x)e Ml ° w x W&-*y V ^ 

and all derivatives in the x direction (i.e., perpendicular to the shock) are considered small 
when compared with terms proportional to k y . This is a perturbation analysis in which the 
transverse wavenumber k y is assumed large and as a consequence the longitudinal wavenum- 
ber ps Wx is also large. The presence of a non-constant amplitude Qx(x) is equivalent to 
keeping the first two terms in an asymptotic expa nsion in the small paramet er (/c y L) _1 . This 



is often called the physical optics approximation (IBender and Orszaglll978l ). 



This analysis is quite standard, but the amusing thing is that we don't even need to 
carry it through. In fact, we shall show later that the stability analysis requires knowledge of 
the physical quantities immediately before the shock, while knowledge of the perturbations 
run with x further from the shock is immaterial. We see from the above that all physical 
quantities close to the shock satisfy 

This is precisely the same form that holds in the homogeneous medium, so that eqs. 1511 and [T71 
still hold. Also, the space-time dependence of the d-mode for the half-plane x > is derived 
in Appendix [Aj assuming spatial homogeneity of the background solution; and this, for the 
argument above, applies in the WKBJ limit also to the upstream fluid. Let us also remark 
that, in this spatially homogeneous limit, vorticity perturbations do not couple to particles: 
in other words, spatial gradients are negligibly small, and thus vorticity perturbations, which 
can only be Lie-advected from upstream infinity, must vanish identically. 

This is the result we need: since we are using a WKBJ approximation, we can treat the 
upstream fluid as if it were homogeneous, with the values for the physical quantities taken 
to be those immediately before the shock: we call this the Homogeneous Approximation: it 
clearly breaks down when the WKBJ analysis does, which occurs for k' 1 ~ L, the size of 
the upstream precursor. 



4. Boundary conditions for the particle distribution function 



Because of diffusion, particles are not restricted to the downstream part of the flow: 
there will be a Sf ^ also upstream, so that we need to consider appropriate conditions to 
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match Sf in the two regions. The first condition is obviously the continuity across the shock, 

Sfi = 8f 2 . (44) 

It is also well-known that the spatial gradient of 5f needs to satisfy a boundary condition 
at the shock: this is derived by integrating eq. H] on an infinitesimal interval straddling the 
shock. The unit vector normal to the surface of the flapping shock is n = (l,tk y Q where ( 
is the shock corrugation amplitude. So, we obtain: 

ft ■ (DVf) 2 - n ■ {DVf\ + h-(v 2 -v 1 ) = 0. (45) 

Writing its first-order linearization and using eq. HH we find: 



D I —la - — li ) + tkyCD 



dx dx 



\dy 2 dy l ) 



+ \ p % {6v2x ~ + {U2 - Ui) = ' (46) 



We note that (df/dy) 2 = (df/dy) 1 = because at zero-th order the medium is uniform in 
coordinates parallel to the shock surface. The result is: 

„{dSf. dSf.\ 1 df , s E , 1 d(5f) . , n 

\dx ~ ~dx J s p i ^ - ^ + aTaT {U2 - Ul) = ° (47) 

Eqs. HH and H7| are the appropriate boundary conditions for our problem. 

We now show how to satisfy the boundary conditions, eqs. HH and H3 at the shock. 
Using the notation introduced in Sect. 12.21 we know that Sf + on the downstream side of the 
shock (the factor e luJt - ik yV will be omitted for simplicity in this subsection) satisfies 

5f+ = A ld 5f ld + g d . (48) 

Please notice that both downstream and upstream the summation is over 2 modes (a pressure 
wave and a third mode). In fact, we know on the one hand that particles perturbations are 
not coupled to entropy and vorticity perturbations, while, on the other hand, we know that 
pressure and third modes, for given values of k y and u>, exist for opposite values of k x (see 
eq. [241 . and thus at least one pressure mode and one third mode will diverge exponentially 
toward infinity: this is unacceptable because we are studying flow instabilities, not variations 
in the boundary conditions at infinity. This discussion will be completed in Section [6j 

On the upstream side we have analogously 



of- — Ai u 5fi u + g u , 



(49) 
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The continuity of 5f at the shock allows us to derive a relationship between the g's: 

9d = 9u + AJfiu - A id 5f id . (50) 



Using 

dSf + 

dx 
d5f. 



—lAidfcxid&fid + ^xd9d — ~ % Aidhxid&fid + ^xd9u + Ai u k x dSfi u — Aidk x dSfid (51) 
1>Ai u k x iu&fiu ~\~ k X u9u (^2) 



<9x 

into eq. H7]we find 

J-){k xd - k xu )g u -\ — = - A iu — -— [pv 2x - 0Vi x ) 

3 amp 3 amp 3 amp 

— DAi U 5fi U (ik X i U + k x d) + D Aid5fid{ik x id + k x d)(53) 
which we regard as an equation for g u (p), whose solution can be written as 

Quip) = Cw c (p) + A iu w iu (p) + A id w id (p) + {Sv 2x - 5u lx )w (p) , (54) 
where the functions w's are derived in the Appendix [B] 

We can now impose the conditions [17] on g u , obtaining: 

C / p 3 vw c dp + A iu / p 3 vw iu dp + A id / p 3 vw id dp + (6u 2x - 5u lx ) / p 3 vw G (p)dp = , (55) 



C J p 3 vw c Ddp + A iu J p 3 vw iu Ddp + A id J p 3 vw id Ddp 

+ (&v 2x - Su lx ) J p 3 vDw (p)($) = , (56) 
And now we can impose the very same conditions on eq. to obtain: 

A id I p 3 v5f id dp - A iu I p 3 v6f iu dp = (57) 



A ld J p 3 v5f id Ddp - A iu J p 3 v5f iu Ddp = (51 



This set of four linear equations in five unknowns (C and the Aj's) can be solved in terms 
of one of them, say A^, the pressure wave downstream. 

We thus see that the conditions at the shock, plus knowledge of the modes, allows us to 
determine the amplitude of all modes (except the vorticity and entropy modes downstream) 
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in terms of the amplitude of the pressure wave downstream. Remembering eq. [3H we now 
see that all fluid quantities at the shock upstream are determined in terms of the amplitude 
of this very same wave; as for the two d-modes, their explicit space- and time-dependence 
is not needed, but is given for completeness' sake in Appendix |A] Now, before discussing 
how to fix u, the shock eigenfrequency, we briefly summarize how to perturb the standard 
Rankine-Hugoniot conditions. 



5. The fluid conditions at the shock 



We follow closely the discussion in lLandau and Lifshita (119871 . Section 90), which re- 



quires only a small adaptation to our problem. In their consideration of the corrugational 
instability, in fact, there were no perturbations on the upstream side of the shock: deviations 
from the unperturbed state were generated by the corrugation of the shock surface (and its 
motion), and led to non-zero perturbations only downstream. In our problem, instead, the 
particles crossing the shock manage to generate new perturbations on the upstream side, 
leading to a slight modification to the Rankine-Hugoniot conditions. From now on, we in- 
dicate with the subscript 1 the quantities on the upstream side of the shock, and with the 



subscript 2 those on the downstream side. Also, to make contact with lLandau and Lifshitz 



work easier, we use as a variable V = 1/p rather than the density p directly. 

We consider a small-amplitude corrugation of the shock surface, away from the x = 
plane, by a small displacement of the form: 

( = ^e**-** (59) 

with respect to which the unit vectors parallel t and normal n to the surface have components 
in the xy plane: 

i=(-lky£,l) S h=(l,lkyQ (60) 

while the surface speed in the direction normal to the surface, with the respect to the 
reference frame of the unperturbed shock, is: 

q-h = iu>( . (61) 

All quantities are, of course, accurate to first order only. 

The first two Rankine-Hugoniot conditions to be perturbed involve the fluid speed, and 
they are the continuity of the fluid speed parallel to the shock surface, and the discontinuity 
of the perpendicular component in terms of perturbed pressure and density (eq. 85.7 of 
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Landau and Lifshitzl 1 19871 ). We have: 



(vi + Sui) ■ i = (v 2 + 



t 



(v t +Sv 1 )-n- (v 2 + Su 2 )-h = y/(P 2 + 5P 2 -Pi- SPJiVi + 8/ 1 -V 2 - SV 2 ) (62) 
whose first-order linearizations are: 

SP 2 - SP X SV 2 - Wx 



Sv 



2,r 



&1„ 



\{V2 -Vx) 



P 2 - Pi Vl - V 2 



(63) 



The next equation to be perturbed is the shock adiabat, which is convenient because it 
is independent of all speeds involved. For a po 

the unperturbed Hugoniot adiabat is given by (ILandau and Lifshitzl 119871 . eq. 89.1) 

V 2 = ( 7 + l)P 1 + ( 7 -l)P 2 
Vx 



whose perturbation yields 

where we have defined 
h 



8V2 
V 2 



47 



ytropic gas like the one w e are considering, 

(64) 

(65) 

(66) 



(7 - 1)P + (7 + 1)P 2 



SVx fSP t 5P 2 



((7 + 1) + (7 - 1)A/Pi)((7 - 1)A/P 2 + (7 + 1)) 
while the ratio P 2 jP\ can be expressed as 

P 2 _ 2 1 M 2 l - (7 - 1) 
Pi 7 + 1 

with Mi the Mach number of the upstream fluid. 



(67) 



We need one more equation: we can use th e equation relating the mas s flux to the 
discontinuities in pressure and density: eq. 85.6 of ILandau and Lifshitzl (119871 ) recites: 

P2-P1 



3 V x -V 2 

where j is the mass flux. When the shock surface is perturbed, the above becomes 



((vi + 



n 



q ■ n) 



P 2 + 5P 2 - Pi - SPx 

(Vx + m) 2 ~ Vi + <5W - V 2 - SV 2 
The perturbation to first order of the above yields 

26u lx 2iuj( 28Vx 5P 2 -5Px SVx - SV 2 



Vl 



Vi 



P 2 - Pi v. - v 2 



(68) 



(69) 



(70) 



Lastly, the condition that 5f be continuous across the shock has already been imposed, 
see eq. HU 
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6. The equation for the perturbation eigenfrequency 

We discuss here first how this problem differs from the one without particles, keeping in 
mind that we are interested in the stability of the shock flapping, so that all modes present 
must have this flapping as their source: we cannot allow modes to come in from spatial 
infinity, because this amounts to perturbation of the boundary conditions, not of the shock 
geometry. 

When we can neglect the particles' pressure, we know that there can be no perturbations 
upstream, since they are either generated at upstream infinity (in which case we would not 
be treating the case of shock instability) or, if generated at the shock, they cannot propagate 
away from it fast enough (the shock is supersonic). In the presence of particles the situation 
changes because particles can diffuse back to the upstream region, so that a perturbation Sf 
generated downstream can return to upstream and perturb the fluid quantities. The situation 
is even more remarkable when one notices that, in this way, there can be a generation of 
pressure waves (even though they are just sonic in a supersonic medium) in the upstream 
region. The reason is shown in eq. [2j the gradient in particles' pressure is a source of 
perturbations, and since particles scatter a finite (i.e., non-zero) distance from the shock, 
there is no obvious reason why even sonic perturbations should not be generated. The 
impossibility of having sonic perturbations in a supersonic medium arises only when the 
point of generation of the perturbations is the shock itself, not a finite distance from it. 

We must also remark that the presence of particles shuffling between downstream and 
upstream, and viceversa, has important consequences also for the kind of waves present 
downstream. In fact, while in the absence of particles the only waves present can be those 
shed by the shock, i.e. entropy, vorticity and pressure disturbances propagating to down- 
stream infinity, when particles are included in the picture we find, by complete analogy with 
the argument above for the upstream region, that they can seed the third mode and the 
d-mode. 

We have seen in Section H] that the amplitude of all modes (and thus of all physical 
quantities) upstream, and of the third mode downstream, can be expressed in terms of a 
single free parameter, the amplitude of pressure waves downstream. Thus all amplitudes 
are fixed, except for three modes downstream (entropy, vorticity, pressure) and the shock 
displacement £; but we also have the four perturbed Rankine-Hugoniot conditions, which 
can be expressed in terms of these four quantities. The vanishing of determinant, once 
substitution of amplitudes for the upstream modes, and for the third mode downstream by 
means of the relations in Section H] has been made, thus fixes the eigenfrequency. 

A word of caution is in order: we have tacitly assumed that the pressure wave propa- 
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gating away from the shock in the downstream region is also the wave with the physically 
acceptable (i.e., decreasing) behaviour at downstream infinity. This is not certain: in fact, 
even in the standard case with no particles, a proper mode exists when the mode explodes 
exponentially in time a nd the pressure wave prop agating away from the shock dies down 
at downstream infinity (ILandau and Lifshitzl 119871 ). This condition can only be checked a 
posteriori, i.e. after having found numerically the value of to. 

We now explicitly derive the set of equations whose determinant must vanish, which is 
what fixes uo for a given k y . 

In the four Rankine-Hugoniot equations there appear the total perturbations of specific 
volume, pressure and velocity. They can be written as sum of the respective perturbation 
for each mode. Downstream we have: 

SV 2 = SV 2s + SV 2p + SV 2t , (71) 
5u 2 = Su 2v + 5u 2p + 5u 2t , (72) 
SP 2 = 5P 2p + 5P 2t , (73) 

where the subscripts s, v, p and t label quantities for entropy, vorticity, pressure waves and 
third mode. Upstream there are neither entropy nor vorticity perturbations: 

SVi = SVip + SVit, (74) 
SJi = Sui p + doxt , (75) 
SP X = 5P lp + 5P lt . (76) 

We must use these expressions in the perturbed Rankine-Hugoniot. Furthermore, we can 
write one of the two components of Su 2v in terms of the other one, through k 2v -Su 2v = 0, where 
k 2vy = k y and k 2vx — uj/v 2 (see eqs. [13]). Then we can simplify our system by eliminating the 
shock displacement (, the remaining component of the velocity perturbation due to vorticity, 
and One equation remains, linking only pressure waves and third modes' perturbations: 

uP 2 (ui - v 2 ) V x {VxPt (uj 2 - k y 2 v lV2 ) 

- [u 2 ((l + h)P 1 -h P 2 ) + ky 2 ((h-l)P x -h P 2 ) v x v 2 ] V 2 ] (5P lp + 8P lt ) 

-cu P x ( Vl - v 2 ) V, {V^ (cu 2 - k y 2 v x v 2 ) 

- [u 2 (hPi -(h-1) P 2 ) + k y 2 (hP x - (1 + h) P 2 ) v x v 2 ] V 2 ) (5P 2p + 5P 2t ) 

+u P x (P x - P 2 ) P 2 (v x - v 2 ) (u 2 - k y 2 v x v 2 ) (V x - V 2 ) (5V lp + 5V lt ) 
-2uPx (Pi - P 2 ) P 2 (to 2 + k y 2 v 2 (-v x + v 2 )) Vx (Vx - V 2 ) (5u lpx + 6u ltx ) 

-2 CO 2 k y Px (Pi - P 2 ) P 2 V 2 Vx (Vx - V 2 ) (SVxpy + Suxty) 

+2 oo 3 Px (Px - P 2 ) P 2 Vx (Vx - V 2 ) (6u 2px + Su 2tx ) 
+2 uj 2 k y Px (Px - P 2 ) P 2 v 2 Vx (Vx - V 2 ) (6u 2py + Su 2ty ) = . (77) 
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Now, eqs. [31] hold both for pressure waves and third modes. We can use them to write the 
above equation in terms of volume perturbations. We remark that each mode has its own z 
and its own k x but, for given k y and uj, they are fixed by the dispersion relation, eq. [2H We 
obtain: 



E TT K(7 - 1) p i + p 2) K - v 2 ) (to 2 - k y 2 Vl v 2 ) (V, - V 2 ) 

4 

l=P,t 

-hj (P 1 -P 2 ) (v 1 -v 2 ) (u 2 + k v 2 v 1 v 2 ) V 2 

+2 (P 1 - P 2 ) (Vi - V 2 ) ((u 2 + fc„ 2 W 2 (V 2 - t>l)) Ziia, + w k y v 2 z liy )\ 

- E tt b p 2 fa - u 2 ) (^ 2 - fc„ 2 vi v 2 ) (v, - v 2 ) 

l=p,t 

-hi (Pi - P 2 ) fa - v 2 ) (to 2 + k y 2 t»i ua) V 2 
+2 (Pi - P 2 ) (Vi - V 2 ) (w z 2ix + k y v 2 z 2iy )} = , (78) 

where the sum is over the pressure mode and the third mode. Now, recalling the definition 
of Ai, we see that A4 = —SVi/V (see also the discussion following eq. En])- Then we have the 
four equations [551 EE EH EH plus eq. [TBI for a total of five linear homogeneous equations in 
five unknowns: four Ai and C. The system has non-trivial solution only if its determinant 
vanishes. This condition determines the eigenfrequency of the system. 



Results 



In this section we ap ply our stability analysis to two shock solutions: a single solution 
(Fig. [H upper panel) by lAmato and Blasil (120051 ) and a multiple solution (Fig. [3], lower 
panel) by lAmato et al\ (120081 ) . 



As we stated above, instability takes place when perturbations exist and grow exponen- 
tially with time. Furthermore, they must decay away from the shock. These conditions may 
be summarized as follows: 



< , %k xld ) < , %k xlu ) > , 
where the subscript i indicates that these conditions must hold for all waves. 



(79) 



For illustrative purposes, first we consider the solution of D'yakov's equation (see lLandau and Lifshitz 
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%J " 2 ( h\ + ^)- (— + ^ 1 (to - u 2 k x ) (l + h L ) = 0, (80) 



19871 . eq. 90.10): 

-)-(- 

u\ V " u 2 ,) \UiU 2 
which fixes the shock eigenfrequency in the linear regime, with 

, 2 dV 2 



h L = r 



dPo 



(81) 

Vi,Pi 



We apply D'yakov's analysis to a test particle solution for a strong shock (upstream 
M ach number Mj — » oo) in a polytropic fluid with index 7 = 5/3. From eqs. 89.6 and 89.9 



Landau and Lifshitzl (119871 ) we obtain, respectively, the compression factor R — V\/V 2 = 4 
and the downstream Mach number l/\/5. The downstream sound speed is c s2 = u 2 /M 2 = 
U\/RM 2 = m 1 v / 5/4. In such system sound waves can propagate away from the shock only 
downstream. The dispersion relation for such perturbations is straightforward: 

L - \ Ul k)j = (k 2 y + kl) , (82) 

where k x is the x-component of the sound wave. To write down the eigenfrequency equation 
for corrugations of the shock we need to calculate h^. This task is straightforward for a 
strong shock because the upstream pressure vanishes. As a consequence, the derivative 
(dV 2 /dP 2 ) Vi p vanishes too (see eq. El]) and hi = 0. Eq. IBOl becomes: 

4cj 3 - \ou\k\ + U\k x [uj 2 + ^ulk^j = . (83) 

Eqs. [82] and E3] form a system to be solved with respect to uj and k x for a given value of k y . 
Since these equations are third-degree homogeneous in lu, k y , k x , if (u,k x ) is a solution for 
a given k y , then (\u>, Xk x ) is a solution for Xk y . Thus the problem is completely solved once 
all the solutions for a given k y are found. Let be k y = 1. Also the upstream fluid speed u\ 
can be set to 1 by redefining the ratio between the units of measure of frequency and wave 
number. We calculated all the solutions of the above equations, for these values of u\ and 
k y , and display them in Tab. [T] 

The first four solutions must be discarded because they correspond to waves propagating 
from downstream infinity to the shock. The fifth solution must be discarded too because it 
diverges exponentially at downstream infinity. More intriguing is the last solution. It seems 
to satisfy all the requirements in order to be a real physical solution and it actually satisfies 
conditions [79] Is it a real instability? No, because it represents a pressure pertu rbation 



advected by the fluid since ui — uk x = 0. From eq. 90.5 in lLandau and Lifshita (119871 ) we see 
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that such a perturbation should have 5P = 0. This is clearly absurd: we assumed k x to be 
the wave vector of a pressure perturbation in order to obtain eq. [80] but we have now found 
a solution with vanishing pressure perturbation. Note that the fifth solution is affected by 
this problem too. 

Summarizing, we solved the equation for the shock eigenfrequency in the test particle 
regime. We found six solutions, but none has physical meaning, since they correspond 
to sound waves either propagating from downstream infinity to the shock surface or with 
vanishing pressure perturbation. 

Below we apply our theory to two particular shock structures (one with multiple solu- 
tion) in order to search for an instability. We do not want to carry out a systematic analysis 
of a set of solutions but we shall just show how our machinery works. For the sake of 
completeness in Fig. [3] we report the particle distributions / we used in our calculations. 

We adopted a system of units of measure so that the following three quantities equal 
1: the speed of light, the proton mass and the numerical constant of the Bohm diffusion 
coefficient, i.e., D p = v(p)p. Below everything will be expressed in these units. 

We proceeded in the following way. We calculated the determinant of the system of five 
eqs. EE [57J ES an d ED] and we set it to 0, obtaining an equation with 5 unknowns: lu, 
k xp u, k xtu , k xpd , kxtd, which are respectively the frequency and the x-components of the wave 
vectors of the upstream pressure wave, the upstream third mode, the downstream pressure 
wave and the downstream third mode. This forms a system of five equations together with 
four dispersions relations as in eq. [22] each one linking u and k y with their respective k x . 
This is exactly what we did above when we solved eqs. [82] and [83] In Fig. H] and in Figs. [5] 
E1I7J we illustrate the solution of our system of equations as a function of k y , respectively for 
the single and the multiple solution of the shock structure. Absolute values of real parts (left 
panels) and imaginary parts (right panels) of to (upper panels), k xpu , k xtu , k xp d, k xt d (lower 
panels) are plotted. Signs of these quantities are reported in Tab. [2] These solutions must 
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Table 1: Solutions of eqs. [821 and [831 for k y = u\ = 1. 
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Fig. 3. — Particle distribution function / as a function of momentum p, measured in units 
of mc. The value of / is divided by the fluid dens ity Vp at upstream infin ity and multiplied 
by p A . Upper panel: single solution, obtained by lAmato and Blasil (120051 ) for a shock with 
Mach number at upstream infinity M\ = 100. Result s of our stabil i ty ana lysis are reported 
in Fig. HI Lower panel: triple solution, obtained by lAmato et all (120081 ) for a shock with 
Mach number at upstream infinity Mi = 100. Results of our stability analysis for particle 
distribution in solid, dotted and dashed lines are reported, respectively, in Figs. El El and [71 
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be all discarded because they have some waves diverging at infinity, just like the first four 
solutions in Tab. HJ We used as starting values in our searching algorithm each one of the 
six corresponding limit solution in the linear case. All the solutions we have found cannot 
be accepted for the same reasons discussed above for the linear regime. 



8. Discussion 

In essence, our method is exact, except for the short- wavelength (WKBJ) approximation 
necessary to treat analytically perturbations in the inhomogeneous u pstream precursor. One 



may however wonder where our method differs from previous work fjMond and O'C. Drury 



19981 ; lToptyginlll999h . which has reported the existence of corrugational instabilities. 



Mond and O'C. Druryl (119981 ) have reported the existence of both genuine corrugational 
instabilities, and of spontaneous emission of sound waves, for some (not all) of their models. 
There is of course a number of differences between this paper and theirs: we do not use the 
two-fluid approximations, and are interested in small-wavelength perturbations, contrary to 
them; we also notice the existence of perturbations in the upstream fluid, which they do not 
discuss. 



Toptyginl (119991 ) included a number of novelties in his treatment, but he too did not 
notice that particles would diffuse upstream, so that he neglects the third and the d-modes 
altogether. Despite this, he does have perturbations upstream (pressure waves!), because he 
remarks that, for sufficiently long- wavelength perturbations, these become supersonic with 
respect to the fluid alone. This occurs because particles and fluid are tightly coupled in 
long wavelength perturbations by diffusion, which traps non-thermal particles, so that the 
restoring pressure is the sum of particles' and fluid's pressures, which is larger than the pure 
fluid speed. This is of course correct, and exists in our computations as well. 

Contrary to these authors, we have found that shocks with non-vanishing particle pres- 
sure are stable to corrugational instabilities, even in the region of parameter space where 
multiple solutions are possible. We believe that the reason for this discrepancy lies in our 
inclusion of diffusion, and the abandonement of the two-fluid approach, as we now argue. We 
stated in the Introduction that the only possibility for the shock destabilization (polytropic 
shocks without particles are well-known to be corrugationally stable) lies in setting up a loop 
whereby perturbations shed by the shock in the downstream region return to the upstream 
region via particle diffusion and excite more perturbations; thus perturbation energy is not 
lost to downstream infinity, but returns to the shock to create more havoc. However, the 
very same mechanism that brings particles back upstream, also causes perturbations' damp- 
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Fig. 4. — Solution of our equations as a function of k y , for the shock structure with accel- 
erated particles' distribution plotted in Fig. [3j upper panel. Upper panels: absolute values 
of real part (left panel) and imaginary part (right panel) of the eigenfrequency to. Lower 
panels: absolute values real parts (left panel) and imaginary parts (right panel) of the x- 
component of the wave vectors k xpu (pressure mode upstream, dotted lines), k xtu (third mode 
upstream, solid lines), k xp d (pressure mode downstream, long-dashed lines), k xt d (third mode 
downstream, short-dashed lines). Signs of these quantity are reported in Tab. [2j 




Fig. 5. — Solution of our equations as a function of k y , for the shock structure with acceler- 
ated particles' distribution plotted in Fig. [3j lower panel, solid line. Please refer to caption 
of Fig. H] for details. 




Fig. 6. — Solution of our equations as a function of k y , for the shock structure with acceler- 
ated particles' distribution plotted in Fig. [3j lower panel, dotted line. Please refer to caption 
of Fig. H] for details. 




Fig. 7. — Solution of our equations as a function of k y , for the shock structure with ac- 
celerated particles' distribution plotted in Fig. [3j lower panel, dashed line. Please refer to 
caption of Fig. @] for details. 
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ening: we have seen that diffusion leads to damping both pressure waves and the so-called 
third modes. Also, we have established that a part of the perturbation occurs in the <i-mode 
(where d stands for damping), where particles are perturbed in phase-space, but the total 
perturbation to the pressure vanishes. Thus, some fraction of the perturbation goes into a 
totally useless form (the G?-mode), the new mode (the third mode) is always damped, and 
even pressure waves acquire a damping which is altogether neglected in the two-fluid approx- 
imation. In the end, while diffusion brings particles (and perturbations) back to the shock, 
the simultaneous damping is so strong to make the excitation of an instability ineffective. 

We should remark that damping of pressure waves is weakest for the largest wavelengths, 
a limit which is unacces sible to us because of o ur W KB J approximation, but is exactly 



the limit investigated by iMond and O'C. Druryl (119981 ). While we have not found a large 



wavelength beyond which the shock becomes unstable, we cannot exclude that a proper 
treatment of the perturbations in the space-dependent precursor may yield a transition to 
the unstable regime. 

A further comment is in order: of all possible geometries, the planar one is probably 
the least likely to display instability. Consider in fact a spherically symmetric explosion like 
a SuperNova. In this case, the perturbations (except of course for entropy and vorticity) 
generated downstream do not escape to infinity, as they do in the planar case, but return 
to the shock because the downstream region is finite and because they are deflected by a 
spatially-dependent refraction index; once they reach the shock, they may generate further 
perturbations. The situation is even more promising when the shock is due to an accretion 
flow, or is stal ling: in fact, except for the presence of particles, this is exactly the scenario 



proposed (see iLamingi 120071 for an analytic approach and discussion) for the generation of 
asymmetries in proto-neutron stars: in this case, the mechanism for the instability of a 
stalled accretion shock is the reflection by the hard star surface into outgoing pressure waves 
of advected entropy perturbations, which return to the shock to generate more mischief. In 
the problem with particles, there is the extra complication due to diffusion, to be o vercome 



to ge nerate instability. Given the relative complexity of even a pure fluid analysis (jLaming 



20071 ). it is likely that this problem will require a numerical approach. 
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A. D-mode in an homogeneous semi-infinite medium 

We give here an explicit expression for the d-mode in an homogeneous but semi-infinite 
medium. The d-mode is the solution of the diffusion equation 

£-V (Al) 

where the speed u is assumed constant because of homogeneity. We solve first the equation 
with u — 0; to do so, we remark that suitable boundary conditions are that 5f — > as x — > 
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±(X), depending on whether we are considering downstream or upstream regions, respectively. 
The solution can be obtained by separation of variables, obtaining: 

Sf = a {p, k x )e vt e k * x e ±lk y y , (A2) 

subject to the constraint 

v = D(k 2 x - hi) . (A3) 

The sign of k x is the one that allows the solution to remain finite at infinity. Solutions 
belonging to different values of k x and p can obviously be superposed, but we know that 
there is another boundary condition (eqs. SHI HH1 and following discussion) to be satisfied: 
at the shock, x — 0, 

Sf = g(p)e mt e- ikyV (A4) 

which obviously gives 

v = iuj , Dk\ = iuj + Dky , a{p) = g(p) . (A5) 
The function g(p) was chosen so that 

p 3 vg(p) dp — D(p)p 3 vg(p) dp = . (A6) 



We have already seen (see the discussion leading to eq. [T7]) that these are the conditions for 
the vanishing of 8P C and its first time derivative at the initial time, and thus at all times. 
The same property is of course acquired by a(p), so that the mode we just found is surely a 
d-mode for the upstream region, in the Homogeneous Approximation. 

When we assume u^O, the above formulae remain correct except for the substitution 
v — > v + uk x . 



B. Explicit derivation of functions w 

We need to find the solution to 

n /, i \ , «2-«i dg u u 2 -u 1 d5f iu 10/ 

D(k xd - k xu )g u -\ — — = - A iu — -— — [du 2x - dv lx ) 

6 amp 6 amp Samp 

+ k xd ) + DA id 5f id (tk xid + k xd )(Bl) 

This is an equation for g u (p)'- 

E(p)g u ( P ) + F^^- = G(p) , (B2) 



E(p) = D(p)(k xd {p) - k xu (p)) , 
_ u 2 - U\ 

F = ""a - ' 

W 3 <91np 3 <91np ; 

-D(p)A iu 8f iu (p)(ik xiu + k xd (p)) + D(p)A id Sf id (p)(ik xid + k xd {p)) 



We find a solution of the homogeneous form of eq. 



d\np 



wc{p) = exp 





= exp 


f- 3 f 


. F Jp m V . 




U 2 ~Ui J p 



rp dri 
(k xd (p')-k xu (p'))D(p')^- 

p 



Now we seek a solution of IB2I of the form g u (p) = C(p)wc(p)'- 

E{p)C{p)w c {p) - ^FC{p)w c {p) + Fw c {p)^^ = G{p) 



dlnp 



Fw c (p)^ = G(p) , 



whose solution is: 



<9 hip 



9u(p) = Cw c {p) + w c {p) 



1 f p G{p') dp' 



Let us define: 



w iu (p) = -w c {p) 



F J Pm w c {p') p' 

p 1 dSf iu (p') dp' 
Pm Wciv 1 ) dlnp' p' 
3 r {%k xiu + k xd (p'))D{p')Sf lu (p') dp' 



w id (p) = w c (p) 
w (p) = -w c {p) 



U 2 - «1 J pm W C (p') p' 

3 r (ik xid + k xd (p'))D(p')Sf ld (p') dp' 



u 2 ~ Mi J Pm w c (p') 
1 f p 1 df(p')dp' 



p' 



u>2 - Mi J Pm w c (p') dYa.pt p' 
therefore we obtain: 

9u(p) = Cw c (p) + A iu w iu (p) + A id w id (p) + (5u 2x - 5u lx )w (p) 
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single solution multiple solution 





Fig. a 


Fig. E 


Fig. [6] 


Fig. [3 
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Table 2: Signs of functions plotted in Figs. EH El [7J None satisfies eq. [791 No instability 
has been found by this analysis. 



